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We study the relation between the eigenfrequencies of the Bogoliubov excitations of Bose-Einstein 
condensates, and the eigenvalues of the Jacobian stability matrix in a variational approach which 
maps the Gross-Pitaevskii equation to a system of equations of motion for the variational parameters. 
We do this for Bose-Einstein condensates with attractive contact interaction in an external trap, 
and for a simple model of a self-trapped Bose-Einstein condensate with attractive 1/r interaction. 
The stationary solutions of the Gross-Pitaevskii equation and Bogoliubov excitations are calculated 
using a finite-difference scheme. The Bogoliubov spectra of the ground and excited state of the self- 
trapped monopolar condensate exhibits a Rydberg-like structure, which can be explained by means 
of a quantum defect theory. On the variational side, we treat the problem using an ansatz of time- 
dependent coupled Gaussians combined with spherical harmonics. We first apply this ansatz to a 
condensate in an external trap without long-range interaction, and calculate the excitation spectrum 
with the help of the time-dependent variational principle. Comparing with the full-numerical results, 
we find a good agreement for the eigenfrequencies of the lowest excitation modes with arbitrary 
angular momenta. The variational method is then applied to calculate the excitations of the self- 
trapped monopolar condensates, and the eigenfrequencies of the excitation modes are compared. 
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I. INTRODUCTION 

In the quantum mechanical description of the ground 
states of Bose-Einstein condensates in the framework of 
the Gross-Pitaevskii equation, the frequencies of elemen- 
tary excitations of the condensates are obtained by solv- 
ing the Bogoliubov-de Gennes equations. In an alter- 
native description, a variational approach with coupled 
Gaussian functions has recently been proposed by Rau et 
al. [U which maps the Gross-Pitaevskii equation to a 
dynamical system for the variational parameters that can 
be analyzed using the familiar tools of classical nonlinear 
dynamics. Ground states correspond to the fixed points 
of the dynamical system, and their stability properties 
follow from the eigenvalues of the Jacobian at the fixed 
points. In this paper we shall investigate the question 
whether or not there is a relation between the eigenvalues 
of the Jacobian and the eigenfrequencies of the quantum 
mechanical Bogoliubov excitations, and if so, to what 
extent the eigenvalues of the Jacobian can reproduce the 
frequencies of these excitations. 

The realization of a Bose-Einstein condensate (BEG) 
with ^^Gr atoms [3] marked the beginning of experimen- 
tal investigations of BEGs with long-range interactions. 
The anisotropic dipole-dipole interaction caused by the 
large magnetic moment of the ^^Gr atoms infiuences the 
properties of the quantum gas and is responsible 
for new phenomena, such as a roton-maxon spectrum 
[5], structured ground states and angular collapse 

[5^. Recently a condensate of ^^''^Dy atoms with an even 
larger magnetic moment was created [9 , lOj, and BEGs of 
other lanthanides with a strong dipole-dipole interaction 
should be possible )1I| . 



A model of a BEG with a different long-range inter- 
action was proposed by O'Dell et al. [12J. In contrast 
to the dipolar interaction, this interaction is monopo- 
lar, i.e., "gravity- like" with an attractive 1/r potential. 
Although it will be difficult to realize this model exper- 
imentally, BEGs with monopolar long-range interaction 
are worth investigating in their own right, since they ex- 
hibit the phenomenon of self-trapping [T^, i.e., the exis- 
tence of a stable condensate without an additional exter- 
nal trap. Furthermore, the isotropic character of the in- 
teraction renders numerical investigations easier than in 
the anisotropic case, and therefore BEGs with monopolar 
interaction can serve as model systems for the treatment 
of condensates with long-range interactions to test new 
approaches and techniques. 

The stationary states of self-trapped monopolar con- 
densates have been calculated in the Thomas-Fermi 
regime and with the variational ansatz of a single Gaus- 
sian |.12j , full- numerically [13j , and with an ansatz of cou- 
pled Gaussians [TJ |2] . Several aspects of the excitation 
spectrum have also been investigated P', [TJl HH], but a 
comprehensive study is still lacking. In this paper we will 
solve the Bogoliubov-de Gennes equations and reveal a 
Rydberg-like structure in the numerically exact Bogoli- 
ubov spectra, similar to the spectra of alkali metals. 

The full-numerical calculations are very accurate, if 
- depending on the method - grid size, number of ba- 
sis functions, etc., are chosen carefully, but may become 
computationally very expensive. As an alternative we 
pursue a variational ansatz with coupled Gaussian func- 
tions. Single Gaussians have been used in the literature 
to obtain qualitative results for BEGs (e.g. in |12| IT6]). 
The ansatz can be extended to time-dependent coupled 
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Gaussians [T71ITS] . and it was demonstrated [IJ[2] that the 
method can quantitatively reproduce the properties of 
the stationary solutions of the Gross-Pitaevskii equation 
with both monopolar and dipolar long-range interaction. 
However, as we discuss below, the ansatz with coupled 
Gaussians can only describe excitations with a maximum 
angular momentum of ^ = 2. Several extensions of a 
Gaussian ansatz have been considered in the literature, 
e.g., Gaussians with Hermite or Laguerre polynomials 
[H Uni [20] , or sines and cosines ^21j . But these methods 
allow for no systematic improvement of the ansatz, which 
is the case for the variational method we present in this 
paper. 

Our variational ansatz is based on a combination of 
coupled Gaussians with spherical harmonics, and can de- 
scribe excitations with arbitrary angular momenta in ra- 
dially symmetric systems. The power of the method will 
be demonstrated by applying it to BECs without and 
with monopolar long-range interaction. 

The paper is organized as follows. In Sec. |TT]we give 
the basic equations, and describe our numerical method 
for calculating the stationary states and excitations of 
self-trapped monopolar condensates. We show that the 
Bogoliubov spectra can be nicely analyzed in terms of 
quantum defect theory. Our variational ansatz with 
time-dependent coupled Gaussians combined with spher- 
ical harmonics is described in Sec. and the equa- 
tions of motion for the Gaussian parameters are derived. 
The method is applied to BECs without and with the 
monopolar long-range interaction. In Sec. |IV| we draw 
conclusions and give an outlook on future work. 



II. FULL-NUMERICAL TREATMENT OF THE 
SELF-TRAPPED MONOPOLAR CONDENSATE 

The time-dependent Gross-Pitaevskii equation (GPE) 
for the self-trapped condensate with short-range contact 
interaction and long-range monopolar interaction reads 
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where a denotes the s-wave scattering length. Since we 
will concentrate on the case of self-trapping, the external 
potential has been omitted. All variables in Eq. ^ are 
given in the natural units introduced in [13j : Lengths 
are measured in units of the "Bohr radius" = /mu, 
energies in units of the "Rydberg energy" — u/2au, 
and time in units of = h/E^. The quantity u is the 
coupling constant of the monopolar interaction defined 
in [12] and depends on the intensity and wave number of 
the laser, and the polarizability of the atoms. 

Eq. ([T|) represents the GPE for the fictitious one-boson 
problem. One can make use of the scaling property of [13 



to scale all quantities to an A'^-boson system: 

(r, a, t, ip) {Nr, N^a, N^t, N-^/^i)). (2) 

The scaled dimensionless units are used throughout this 
work and in all figures whenever considering monopo- 
lar condensates. In these units, the only remaining pa- 
rameter is the scattering length a The station- 
ary GPE can be obtained by substituting ilj{r,t) = 
ip{r) exp(—ifit), with the chemical potential in the 
time-dependent GPE ([ij, which leads to 
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A. Calculation of stationary solutions 

For a numerical treatment of the stationary GPE ^ 
it is convenient to convert the integro-differential equa- 
tion into two coupled differential equations. This can be 
achieved by defining the mean-field potential 



0(r) 
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Since we search for radially symmetric stationary so- 
lutions we assume the wave function and the mean- 
field potential to depend only on the radial coordinate: 
■0(r) = ip^r) and 0(r) = 4){r). Letting the Laplacian in 
spherical coordinates act on Eq. Q one obtains the two 
one-dimensional, nonlinear coupled differential equations 
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(r)-87r|V'(r)|^ =0. (5b) 



The system of Eqs. ^ can be solved numerically in dif- 
ferent ways. Since it is a one-dimensional problem, one 
can integrate the equations using a Runge-Kutta algo- 
rithm from r = to a sufficiently large value rmax with 
appropriately chosen initial conditions for ip{0), V''(0), 
(/)(0) and (j}'{0) [H [TSl US]. Their values must be var- 
ied until the wave function converges towards zero at 
f — ''max- With this method the ground and excited 
state can be calculated efficiently. However, to obtain 
a normalized solution ip{r) the wave function, scattering 
length, and mean field energy must be rescaled. Thus, it 
is difficult to obtain a solution for a given fixed value of 
the scattering length. Additionally, it is not easy to cal- 
culate the Bogoliubov spectrum of the system with this 
method, since the solutions of the Bogoliubov-de Gennes 
(BDG) equations have large extensions, and a very big 
value of Tinax has to be chosen. For example, to calcu- 
late 20 eigenvalues for an angular momentum of / = 6, 
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Tmax needs to be larger than 1000. In this case, machine 
precision in the Runge-Kutta method is not sufficient to 
obtain converged solutions, leaving this method useless 
for higher modes. In l^, TF, only the three lowest I ~ 
modes could be calculated. 

Another method is the imaginary time evolution (re- 
placement t t = IT in Eq. Q) of an initial wave func- 
tion on a grid. As time evolves the wave function con- 
verges to the ground state. This method is useful to find 
the ground state or a metastable state of a system. How- 
ever, a collectively excited state, as we consider below, 
cannot be obtained by imaginary time evolution. 

To avoid these disadvantages, we use the finite- 
difference method to solve the stationary GPE (|5| : Wave 
functions and the mean-field potential are discretized on 
a grid and all derivatives are replaced by their finite- 
difference approximation. To arrive at a closed system of 
algebraic equations which can be solved by a nonlinear 
root search one needs appropriate boundary conditions: 
ip'{0) — and (p'iO) — 0, to ensure that the functions 
are differentiable at the origin, and V'(''max) = to ob- 
tain a normalizable wave function. The fourth boundary 
condition can be obtained by looking at the asymptotic 
behavior of the mean-field potential Q . Approximating 
1/ |r — r'l « 1/r for r ^ r' and assuming a normalized 
wave function ^, one obtains from Eq. Q 0(r) ss — 2/r. 
The fourth boundary condition is therefore set to be 

We perform the nonlinear root search using the Powell 
hybrid method |23| . In addition to the equations origi- 



nating from the finite-difference scheme, the normaliza- 
tion condition has to be included, as well as the chemical 
potential as a parameter to be determined by the root 
search. 



B. Bogoliubov-de Gennes equations 

The stability and elementary excitations of a self- 
trapped monopolar condensate have already been an- 
alyzed in the literature: the lowest monopole and 
quadrupole oscillation analytically and numerically |14| . 
the first monopole modes [52], and the lowest monopole 
and quadrupole modes by means of a variational ansatz 
with coupled Gaussians [2 . However, to the best of our 
knowledge, a calculation of the Bogoliubov spectrum for 
arbitrary angular momenta and higher excitations does 
not yet exist. 

To derive the BDG equations, one starts from the usual 
ansatz for a perturbation of a stationary state 
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where cj is the frequency and A the amplitude of the per- 
turbation (|A| <C 1), and /i is the chemical potential of 
the stationary solution with corresponding mean-field 
potential 0o- Eq- ^ is inserted into the time-dependent 
GPE ([ij, terms of second order in A are neglected, and 
collecting terms evolving in time with exp(— iwi) and 
exp{iu!t) yields the BDG equations 
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with the auxiliary field (similar to the mean-field poten- 
tial) 



/(r) = -2j dV' 
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The ansatz of Eq. (l6| possesses a symmetry: the ex- 
change of u{r) -H> v*yr) and uj o —uj leaves the ansatz 
invariant. Thus for each solution (u, v) and a; of Eqs. ([7|, 
{v* ,u*) with —UJ is another solution and both solutions 
represent the same physical motion. For that reason, only 
solutions with Rew > need to be considered. There are 
two solutions of Eqs. ([t]) which deserve special attention. 
If "00 is assumed to be real, then u{r) ~ ~v{r) — ipo{r) 
is a solution of Eqs. (|7| with the frequency uj — 0. This 
represents the well-known gauge transformation of the 
condensate wave function 'ijj{r) — > ^/;(r) exp(i0) with a 
real phase 0. This gauge mode does not describe a phys- 



ical motion of the condensate, and since it is always part 
of the Bogoliubov spectrum, we will not discuss it when 
presenting the results. 

Furthermore, there always exist solutions of the BDG 
equations with frequencies identical to the trapping fre- 
quencies [2^. These modes represent the center-of-mass 
oscillations of the condensate along the three space direc- 
tions with angular momentum 1 = 1. In the case of the 
self-trapped monopolar condensate, there are no external 
traps and therefore the frequencies are w = 0, which cor- 
responds to a constant displacement of the condensate. 

Since the wave function "00 £^nd the mean-field poten- 
tial 00 s-re radially symmetric, we can separate the solu- 
tions u and V by means of spherical harmonics 
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with the radial (excitation) quantum number n and the 
usual angular momentum quantum numbers l,m. Us- 
ing the multipole expansion of the integration kernel 
l/\r — r'\ (see, e.g., [2F and Eq. (A.29l), we can also 
express the auxiliary field ^ in the form fnimi'r) — 
Yim{(^,<P)fni{T) with {tpQ and 0o £^re assumed to be real 
from now on) 



where r< = min(r, r') and r> = max(r, r'), respectively. 
Inserting the Laplacian in spherical coordinates and using 
the separation ([9]), we finally obtain from Eqs. Q 
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We solve Eqs. (Ill using the finite-difference method. 



After choosing a grid, approximating the derivatives by 
finite differences and replacing the integral in the auxil- 
iary field / by an appropriate integration rule (we use the 
trapezoidal rule), Eqs. ( 11 1 turn into a matrix eigenvalue 
problem 



M 
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The eigenvalues of the matrix M can then be found by 
numerical diagonalization. 

In actual calculations we found it advantageous to 
choose a non-equidistant grid, since the solutions u and 
V can be highly oscillatory near the origin, and at the 
same time extend to large values of r. We use partially 
equidistant grids, i.e., an equidistant grid with step size 
Ari in the interval [0,ri], another equidistant grid with 
a different Ar2 in the interval [ri,r2], etc. 



C. Results 

Since the properties of the stationary solution have 
been discussed in detail in the literature |13l [T51 E^ . we 
only give a brief review. Our results coincide with those 
obtained using the outward integration method, and thus 
for the stationary states both methods can be considered 
equally applicable. In Fig. [TJwe plot the mean-field en- 
ergy Erai and the chemical potential ^ of the ground and 
excited state as a function of the scattering length a. Two 
solutions are born in a tangent bifurcation at the critical 
scattering length a = Ocrit ~ —1.025. At this point, the 
mean-field energy, chemical potential and wave functions 
of the ground and excited state merge. For a — 0, the 
mean-field energy and chemical potential of the excited 



00 



4 ^ . . . 

ground state 

3 ' excited state 

2 - 
1 - 



_i L, 1 1 1 — 

I 

-1 ■^'^^^^^^^^ZiZ 

-2 

-3 - 

-4 - 

-5 - 

-6 - ground state 

~7 ' excited state 

-8 ^ ' ' ' — 

-1 -0.8 -0.6 -0.4 

scattering length a 



-0.2 



Figure 1. (Color online) Mean-field energy E^^i and chemical 
potential /i of the ground and excited state of a self-trapped 
monopolar condensate as functions of the scattering length 
a. For a scattering length lower than the critical value of 
aciit ~ —1.025 no stationary solution exists. At a = Ocrit 
the two solutions emerge in a tangent bifurcation. For the 
ground state, both Emi and fi stay negative in the range of 
the scattering length considered. These quantities diverge for 
the excited state in the limit o — >■ 0. 



state diverge, implying that this state does not exist for 
a > 0. 

Using the method described in Sec. |IIB| we have calcu- 
lated the Bogoliubov spectrum of the ground state. For 
the angular momenta from ^ = to 3, Fig. [2] shows the 
frequencies of the Bogoliubov excitations as a function 
of the scattering length a. The ground state is stable, 
since its spectrum contains only real frequencies. It can 
be seen that as the scattering length is decreased towards 
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Figure 2. (Color online) Frequencies of Bogoliubov excitations 
of the ground state in Fig. [l] for the angular momenta from 
i = to 3 as functions of the scattering length a. The seven 
lowest eigenvalues are shown for each angular momentum. 
The spectrum only contains real frequencies, i.e., the ground 
state is stable. The lowest mode for I — tends to zero as 
a — >■ Ocrit which leads to the collapse of the condensate. The 
lowest I = 1 mode corresponds to a displacement of the center- 
of-mass of the condensate, while its shape remains unaffected. 
The frequency of this mode is exactly the trapping frequency 
|24| . in this case lj = 0. The frequencies of the other modes 
increase as the scattering length is decreased, finally merging 
with the modes of the excited state for a — ^ flcrit (see Fig.[3|. 
Note that for fixed scattering length the distance between 
adjacent frequencies diminishes with growing radial quantum 
number, indicating the convergence of the frequencies to a 
(scattering length dependent) limit frequency. 



its critical value the frequency of the lowest mode with 
Z = at first slightly increases but then goes to zero at 
a — >■ Ocrit, where the state vanishes. This mode is re- 
sponsible for the collapse of the condensate. The lowest 
I = 1 mode has the frequency uj — and corresponds 
to the displacement of the center-of-mass of the conden- 
sate. This frequency remains constantly lu — as the 
scattering length is varied, and, as already mentioned, 
corresponds to the (vanishing) trapping frequency. 

The results for the excited state are presented in Fig.|3] 
All frequencies merge with those of the ground state 
modes at the critical scattering length. There exists one 
imaginary frequency for the angular momentum 1 = 0. 
Therefore the excited state is unstable with respect to 
this excitation, which leads to a collapse with / = sym- 
metry. As for the ground state the lowest mode with 
I = 1 represents the displacement of the condensate and 
is constantly oj = 0. 

In Fig. |4]the Bogoliubov functions u and v are shown 
for the angular momentum 1 = 0. The lowest functions 
with n = 1 and n = 2 are concentrated near the ori- 
gin and have the same extension as the wave function 
of the stationary solution (see Fig. |6]). For the higher 
modes, the functions u extend further out, which is a 
consequence of the missing external trapping potential. 
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Figure 3. (Color online) Same as Fig. [2j but for the excited 
state. There is one imaginary frequency for 1 = 0: the ex- 
cited state is unstable with respect to small perturbations. 
As for the ground state, the lowest I = 1 mode is a; = 
and corresponds to a displacement of the center-of-mass of 
the condensate. Again, the frequencies of the stable modes 
apparently converge to a limit for fixed scattering length. 
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Figure 4. (Color online) Bogoliubov functions u„i{r) and 
Vni{r) for the angular momentum I = and the radial quan- 
tum numbers n = 1, . . . , 5. The scattering length is a = —0.4. 
The mode with n = 1 represents the gauge mode discussed 
in Sec. II B and the functions Uio(r) and Vio{r) are equal to 



the stationary solution ipo, except for the sign. The functions 
Unoir) have n — 1 nodes, whereas all functions v„o{r) show 
qualitatively the same behavior and are nodeless for all n. It 
can be seen that with growing radial quantum number the 
functions «„o(f) extend to ever increasing values of r. 



D. Quantum defect analysis of the Bogoliubov 
spectrum 

To prove that for given scattering length the frequen- 
cies of the Bogoliubov excitations converge to a limiting 
frequency we determined the 20 lowest modes for the 
angular momenta Z = to 6. As an example. Fig. [5] 
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Figure 5. (Color online) Frequencies of the Bogoliubov excita- 
tions of the ground state of a self-trapped monopolar BEC for 
a fixed scattering length a = —0.4, plotted for different values 
of the angular momentum. The dotted line gives the value of 
the chemical potential. Obviously, as observed in Fig. [2] and 
Fig. [3] the frequencies converge to a common limit, which 
is the chemical potential. The Bogoliubov spectrum can be 
described by a Rydberg formula with quantum defects. 
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Figure 6. (Color online) (a) Wave function ipo and (b) 
mean-field potential 00 for the ground state of a self-trapped 
monopolar condensate at a scattering length of a = —0.4 as 
functions of the radial coordinate r. The wave function ap- 
proaches zero exponentially, whereas the mean-field potential 
behaves like — 2/r for large values of r. In this region, the 
wave function can be neglected and the mean-field potential 
replaced by its asymptotic form in the BDG equations. 



shows, for the scattering length a = —0.4, the Bogoli- 
ubov spectrum of the ground state. The convergence of 
the frequencies to a common limit, independent of Z, is 
evident. The spectrum is reminiscent of Rydberg spectra 
known from alkali atoms. Similar to the spectra of these 
atoms, the structure of the Bogoliubov spectra can be 
understood in terms of quantum defect theory. 



plify due to the fact that the wave function decays ex- 
ponentially, and the mean-field potential converges to 
(j)o{r) fa —2/r (see Fig. [6]). Setting '(/'o(^) ~ for r > r^, 
all terms containing ipo can be neglected in (11 1, and (/jq 
can be approximated by —2/r. This leads to the asymp- 
totic form of the BDG equations 



-UJnlVnlir) = 



dr^ 



2d l{l + 1) 
r dr 



d^ 2d l{l + 1) 
dr^ r dr r^ 



Uniir), 
(13a) 

v„i{r). 
(13b) 



Obviously in this limit u and v obey the same equa- 
tion, namely the Schrodinger equation of the Coulomb 
problem, except for the opposite sign of the eigenvalues. 
Therefore asymptotically only one equation of ( 13 ) needs 
to be considered (which will be the one for u) . The scat- 
tering length enters into Eqs. (13 1 only indirectly via 
fi = ii{a). 

The approximations made are only valid, if the func- 
tion values of u and v are small for r < rc. Especially 
for lower angular momenta this is not the case. In the 
physics of alkali metals a similar problem occurs: The 
valence electron far away from the nucleus "feels" an at- 
tractive — 1/r potential, which results from the shielding 
of the core electrons. Close to the nucleus, the core elec- 
trons and the true nuclear potential has to be considered. 
A similar situation happens here, cf. Fig. [6] To account 
for the deviation of the potential from the pure Coulomb 
potential at smaller values of the radial coordinate we 
can also introduce a quantum defect in the formula for 
the Rydberg series eigenvalues (see, e.g., |26|). 
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w„i = -11 - 



(14) 



where the quantum defects 5i depend on the angular mo- 
mentum. The negative chemical potential is the limit of 
the frequencies for n — > oo. The quantum defects can be 
obtained by least-squares fits of the Bogoliubov frequen- 
cies ujni to Eq. (14). They converge to constant values 
for large n. Since Eq. ( 14 1 strictly holds only in this 



For large values of r the BDG equations (11) sim- 



limit, in the fits it can be necessary to neglect the lowest 
frequencies. 

For growing angular momentum, the repulsive effec- 
tive potential l{l + l)/r^ becomes stronger, and this cen- 
trifugal barrier ensures that the absolute values of the 
functions u and v decrease close to the origin r = 0. 
This leads to a smaller quantum defect (5;, since the ap- 
proximation made in deriving Eqs. (13) becomes valid 
at smaller values of r. In accordance with the quantum 
defects in alkalis the quantum defects therefore will 
tend to zero for higher angular momenta. 

In Fig. [7] we present the quantum defects calculated 
for the Bogoliubov excitations of the ground state. Ob- 
viously the quantum defects for I — and I = 1 show 
a strong dependence on the scattering length, while for 
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Figure 7. (Color online) Calculated quantum defects 5i for 
the Bogoliubov spectrum of the ground state (Fig. [2| in de- 
pendence of the scattering length a for different angular mo- 
menta I. The quantum defects for / — and / — 1 rise steeply 
and turn from negative to positive as the scattering length 
is decreased, while the quantum defect for / = 2 shows only 
a weak dependence on the scattering length and drops close 
to the critical scattering length. As expected, for the higher 
angular momenta / > 2 the quantum defects are close to zero. 



I > 2 they are almost constant, and in particular close to 
zero for I > 2. Eq. (14 1 reproduces the frequencies of the 



Bogoliubov excitations of the ground state for all modes 
with an absolute error of less than 10^^, except for the 
two lowest I = modes and the lowest I — 1 mode. The 
quantum defect analysis for the Bogoliubov excitations of 
the excited state is presented in Fig. |8] The quantitative 
statements made for the excitations of the ground state 
also hold for this state. The only difference is that the 
quantum defect for I = 2 tends to zero as the scattering 
length is increased. 

Thus by means of quantum defect analysis we have 
been able to explain the Rydberg-like structure of the 
Bogoliubov spectra of the ground and excited state of 
self-trapped monopolar BECs, and could confirm that 
the negative chemical potential is the limit of the fre- 
quencies for all angular momenta. 



III. VARIATIONAL APPROACH WITH 
GAUSSIAN FUNCTIONS AND SPHERICAL 
HARMONICS 

We now turn our attention to variational calculations. 
The simplest ansatz with a single Gaussian centered at 
the origin was used by Perez-Garcia et al. fTF| to deter- 
mine monopolar and quadrupolar modes of BECs with- 
out long-range interactions. The ansatz was improved by 
using coupled Gaussians [171 15] , and it was shown [T] |2] 
that this method is capable of reproducing accurately the 
stationary states even of BECs with long-range interac- 
tions, calculated numerically. The ansatz employed to 



Figure 8. (Color online) Same as Fig. [Tj but for the excited 
state. Since the Bogoliubov spectra of the ground and excited 
state merge at the critical scattering length, the same holds 
for the quantum defects. The quantum defects So and 5i 
grow as the scattering length is increased, whereas S2 drops 
and tends to zero. As in the case of the ground state, the 
quantum defects for I > 2 are close to zero. 



determine the stationary solution of a radially symmet- 
ric condensate was 



N 



(15) 



fc=i 



where the complex quantities Aj! and 7'"' are the widths 
and the amplitudes, respectively, of each Gaussian. The 
above ansatz can only describe monopolar excitation 
modes, since the wave function ip is independent of the 
angular coordinates 9 and (j>. If one chooses the widths 
differently for each space direction. 



N 

fc=i 



(16) 



the width of a condensate can oscillate independently in 
each direction, which represents quadrupolar oscillations. 
A generalization of Eqs. (15 1 and (16 1, which includes 



general square and linear terms in the exponentials, is 

HZIITH] 



N N 
k=l k=l 



= ^ exp [-r^A'^r - {p'')^r - ■ 



(17) 



with complex symmetric matrices A'^, complex vectors 
p'' and complex numbers 7^^. This ansatz can describe 
excitation modes with angular momenta up to ^ = 2. To 
see this consider a small deviation Sz of the variational 
parameters from those of a stationary solution Zq and 
Taylor expand the ansatz with coupled Gaussians (17 1 
for the perturbed wave function ipi^Q +^-2) to first order 



8 



the matrix K and the vector h are defined by 



dip = 6z ■ 



dip 
dz 



N 

= -J2 (r^<5AV + {Sp'^fr + dj'') g\^^^ . (18) 
fc=i 

Since only terms at most quadratic in x, y, z appear in 
front of the exponentials, these terms can be expressed 
by spherical harmonics with angular momenta Z = 0, 1, 2, 
which proves our statement. 

We apply an ansatz which is capable of describing exci- 
tations with - in principle - arbitrary angular momenta. 
Motivated by the separation in the BDG equations with 
spherical harmonics in Eq. ([9]), we directly include the 
spherical harmonics in an extended variational ansatz 



N 



k=l 



(19) 



The amplitudes account for additional angular mo- 
menta {l,m). The sum over {l,m) may include arbitrary 
angular momenta, adjusted to the problem. For instance, 
if one wishes to calculate the linear perturbation of a spe- 
cific angular momentum {l,m), as we do below, the sum 
in Eq. ( 19 1 needs to include [l, m) and {I, —m), since the 
nonlinear terms in the GPE lead to a coupling of different 
angular momenta. 



A. Equations of motion and stability analysis 

In order to carry out calculations with the extended 
variational ansatz ( 19 1 , we need the equations of motion 



for the variational parameters. We use the approach of [1' 
based on the Dirac-Frankel-McLachlan time-dependent 
variational principle [271 [28] . An arbitrary ansatz for the 
wave function is made ip = ip{z), with the - in general 
complex - variational parameters z — (zi, . . . , zm), for a 
system governed by the Schrodinger equation 



itp — Hip, 



(20) 



where the Hamiltonian H may contain nonlinear terms 
in the wave function. The principle states that the norm 
of the difference between the left- and the right-hand side 
of (pO| 



i^\m)-Hijj{t)\f 



(21) 



must be minimized. For a fixed time t, ip{t) is given, 
and / is minimized by varying (p{t). After the minimiza- 
tion, (p is set to (p = ip. A necessary condition for the 
minimization of / is |15] 



Ki 



-i/i, 



(22) 



dip 


dip 


dzi 


dzj 



hi = 



dip 
dzi 



Hip 



(23a) 
(23b) 



Stationary solutions can then be found by requiring 



M 



(K-i) 




(24) 



and searching for z in a, nonlinear root search. 

The stability properties and linear oscillations of a 
stationary solution Zq can be found by first changing 
from the complex M-dimensional vector 2; to a real 2M- 
dimensional vector z containing the real and imaginary 
parts of the variational parameters, and considering a 
small perturbation, z{t) = Zq + Sz{t). Linearization of 
the equations of motion ( p2| yields the time dependency 
of the perturbation [T] 



with the Jacobian 



Si{t) = JSz{t) 



dzj 



(25) 



(26) 



evaluated at the fixed point i = io- The excitation 
modes of the stationary solutions are finally found by 
diagonalizing the Jacobian J. 



All integrals appearing in Eq. ( 22 1 with the ansatz ( 19 1 
can be calculated analytically. The contact interaction 
leads to integrals over four spherical harmonics which 
can be expressed in terms of Wigner-3j symbols. The 
contribution of long-range monopolar potential can be 
evaluated by inserting the multipole expansion for the 
monopolar integration kernel, which leads to Gaussian 
integrals. For further details of the calculations we refer 
to the appendix. 



B. Test in a system without long-range interactions 



As a first test we apply the extended variational ansatz 



( 19 1 to a condensate in a radially symmetric trap with 
short-range scattering interaction. The GPE for this sys- 
tem reads 



■ dip , 



-A + r^ + 8Tra\iP{r,t)f ip{r,t). (27) 



Here units based on the trapping frequency "f — uj/2 and 
the harmonic oscillator length oq = ^Jujnvjj have been 
used. The scaled dimensionless scattering length a in 



(27 1 corresponds to Na/a^) in SI units, with the particle 
number N . These units will be used in all figures for the 
condensate without long-range interaction. The BDG 
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Figure 9. (Color online) Comparison of the full-numerical Bo- 
goliubov spectrum of a BEC with attractive contact interac- 
tion with the spectrum obtained from the variational ansatz 
with coupled Gaussians and spherical harmonics (SH). The 
variational ansatz has been used with 5 coupled Gaussians 
and spherical harmonics up to an angular momentum of Z = 3. 
For the lowest modes we find excellent agreement. There are 
almost no deviations for frequencies u < 10. Just slightly 
above the critical scattering length small differences can be 
seen in the figure. For the higher modes, differences become 
larger and the variational ansatz can describe the Bogoliubov 
modes only qualitatively correct. 



equations are given in Eqs. ([?]) and ( [TT| , respectively, 
with all terms containing the mean-field potential (/jq and 
the auxiliary field / omitted, and the trapping potential 
Vext = included. 

The BDG equations for condensates with short-range 
interaction were first solved numerically by |29| [30] . In 
this work, we used the method discussed in Sec. |II B[ 

Fig. [9] shows the eigenfrequencies of the Bogoliubov 
excitations of the ground state with I ~ 0,1,2 and 3 
as functions of the scattering length. For a — one 
obtains the equidistant eigenfrequencies of the harmonic 
oscillator. When the scattering length is decreased the 
attractive short-range interaction acts as a perturbation, 
and the frequencies are slightly shifted. For a — >■ Ocrit ~ 
—0.0575 the lowest I = mode drops to zero marking the 
collapse of the condensate. The lowest mode with / = 
1 represents the oscillation of the center-of-mass of the 
condensate, and its value is exactly that of the trapping 
frequency a; = 27 = 2. 

For comparison in Fig. |9]we also show the results for 
the eigenvalues of the Jacobian matrix at the ground 
state fixed point obtained in the variational ansatz us- 
ing 5 Gaussians in combination with spherical harmonics 
(19 1. One recognizes that in particular the eigenvalues 



of the lowest modes in each angular momentum band 
excellently agree with the eigenfrequencies of the Bogoli- 
ubov excitations. It is only close to the critical scattering 
length that small deviations appear. The lowest center- 
of-mass excitation with I = 1 can even be reproduced 



Figure 10. (Color online) Comparison of both spectra as in 
Fig. |9j but here for a fixed scattering length of a = —0.4 
and angular momenta up to / = 6. For / = the variational 
ansatz reproduces the Bogoliubov frequencies very well for 
the four lowest modes, and with only small deviations for the 
two lowest modes in the higher angular momentum bands. 



within numerical accuracy, independent of the number 
of Gaussians used. For the higher modes with eigenval- 
ues of the Jacobian lu > 10, only far away from the crit- 
ical point the variational and full-numerical results still 
approximately correspond to each other, and in the vicin- 
ity of the critical scattering length the Jacobi eigenvalues 
can reproduce the behavior of the Bogoliubov excitation 
eigenfrequencies only qualitatively. 



We also tested the variational ansatz ( 19 1 for higher 
angular momenta up to Z = 6. The results for a fixed scat- 
tering length of a = —0.4 are presented in Fig. [TO] One 
recognizes a very good agreement for the lowest modes in 
each I band, and small differences for the second-lowest 
modes. This demonstrates that for condensates with at- 
tractive short-range interaction the eigenvalues of the Ja- 
cobian matrix calculated at the fixed point corresponding 
to the ground state in the new variational ansatz indeed 
quantitatively coincide with the eigenfrequencies of the 
lowest Bogoliubov modes. 



Application of the variational approach to the 
monopolar condensate 



We now apply the extended variational ansatz (19 1 to 
the self-trapped monopolar condensate. For the three 
lowest excitations Fig. [TT] shows the comparison of the 
full-numerical Bogoliubov spectrum with the spectrum 
obtained from the eigenvalues of the Jacobian matrix in 
the variational ansatz. We used = 6 Gaussians and 
angular momenta up to ^ = 3. The lowest modes for I — 
and I — 1 match very well in the whole range of scatter- 
ing lengths considered. For the second-lowest I — and 
I = 1 and the lowest / = 2 mode we find a good agree- 
ment, but the differences become larger as the scattering 
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Figure 11. (Color online) Comparison of the full- numerical Bogoliubov spectrum of the ground state of a self-trapped monopolar 
BEC with the spectrum obtained from the variational ansatz with coupled Gaussians and spherical harmonics (SH). The 
variational ansatz has been used with 6 coupled Gaussians and spherical harmonics up to an angular momentum of / = 3. For 
the lowest I = and I — 1 mode the results of both methods almost cannot be distinguished. The differences of the frequencies 
of the second lowest I — and I — 1 and the lowest / = 2 modes are small in the range of the scattering length considered. The 
lowest / = 3 mode is well approximated by the variational ansatz, but the differences in frequency are larger compared to the 
frequencies belonging to lower angular momenta. For the higher modes there is no quantitative agreement. 



length approaches the critical point. Nevertheless, we 
have the result that the variational ansatz with coupled 
Gaussians and spherical harmonics is a valid alternative 
to the full-numerical quantum mechanical approach also 
in this case, if one is interested in these modes. 

Looking at the lowest mode with / = 3 one finds that 
the agreement is good for scattering lengths around a = 
0, but the two frequencies deviate as the scattering length 
is decreased. The eigenmode of the variational ansatz can 
only be seen as an approximation to the full-numerical 
one. The other modes can only be described qualitatively 
by the variational approach. 

We also applied the variational ansatz p9| for higher 
angular momenta up to Z = 6. The results for a fixed 
scattering length of a = —0.4 are presented in Fig. [12] 
As already noticed, only the lowest modes and angular 
momenta agree well with the numerically exact values. 
In the remaining modes the excitation frequencies differ 
distinctly. For I = 5, the frequency of the lowest mode 
even lies above the negative chemical potential, which 
could be identified as the upper limit of the Bogoliubov 



spectrum. Obviously, the variational ansatz with cou- 
pled Gaussians and spherical harmonics is not as appro- 
priate for the self-trapped monopolar condensate as for 
the condensate without long-range interaction. To ob- 
tain still better results in the variational ansatz, it would 
be necessary to use more than = 6 coupled Gaussians. 
This, however, leads to numerical difficulties, since the 



matrix K in the equations of motion ( 22 I becomes more 



and more ill-conditioned when the number of Gaussians 
is increased, which leads to an inaccurate solution of the 
linear system of equations. 



IV. CONCLUSION AND OUTLOOK 

We investigated the Bogoliubov spectrum of the self- 
trapped monopolar condensate full-numerically with the 
finite-difference method. With this method, we were able 
to calculate many modes for angular momenta from I ~ 
to / = 6. We found a similar structure as in the spectra of 
alkali atoms. The behavior could be explained by quan- 
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Figure 12. (Color online) Comparison of both spectra as in 
Fig. |11| for a self-trapped monopolar BEC, at the fixed scat- 
tering length of a = —0.4 and angular momenta up to Z = 6. 
For I = and I — 1 the two lowest modes, and for I = 2 and 
I — 3 only the lowest modes, agree well. For higher angular 
momenta I > 5, the lowest mode lies even above the limit of 
the numerical Bogoliubov spectrum (compare with Fig. [sf . 



turn defect theory, and it was found that practically the 
entire spectrum can be described by a simple Rydberg 
formula with quantum defects. 

As an alternative to full-numerical calculations of con- 
densate excitations a new variational ansatz was pre- 
sented which combines coupled Gaussians with spherical 
harmonics. Using the time-dependent variational princi- 
ple we could derive the equations of motion for the vari- 
ational parameters. We applied the variational ansatz 
to two different systems. For condensates with an at- 
tractive short-range interaction we found that there is a 
good agreement between the quantum mechanical eigen- 
frequencies of the lowest Bogoliubov excitations and the 
eigenvalues of the Jacobian stability matrix. In this way 
we have been able to link the concepts of stability in 
quantum mechanics and in classical dynamical systems 
with each other. 

For self-trapped condensates with additional 1/r in- 
teraction we also found a good agreement for the very 
lowest modes, but the variational ansatz works less well 
for higher modes. What is the reason for this? For the 
condensate without long-range interaction in an external 
trap, the confining radially symmetric harmonic potential 
dominates the properties of the system in a wide range 
of the scattering length. The contact interaction quasi 
acts as a perturbation. Therefore, a variational ansatz 
in which the radial part is determined by Gaussians is 
very well adapted to describe the stationary solutions 
and their excitations. 

For the self-trapped monopolar condensate, on the 
other hand, an external trap is missing and the inter- 
actions alone determine the properties of the system. As 
pointed out in Sec. |IID| the asymptotic form for r — >■ oo 
of the BDG equations is equivalent to the Schrodinger 



equation of the hydrogen atom. Therefore in that range 
the solutions u and v could be approximated by Laguerre 
polynomials and the exponential function exp(— ar) with 
some a > 0. A variational ansatz with coupled Gaussians 
and spherical harmonics is not well suited to reproduce 
this asymptotic behavior. However, as soon as a radially 
symmetric trap is switched on, the agreement between 
the quantum mechanical and the nonlinear dynamics ex- 
citations is present again also for the higher modes. 

All together it was shown that especially in the case 
without long-range interactions the extended variational 
ansatz works well and can reproduce the lowest modes for 
arbitrary angular momenta, which is significant progress 
compared to the ansatz with coupled Gaussians only. If 
one is interested only in the lowest modes, the ansatz is 
a valid alternative to the full-numerical calculations. 

So far, we only calculated the linear dynamics in the 
vicinity of a stationary solution. It remains to be shown 
whether or not the ansatz is capable of describing also 
the full nonlinear dynamics of a BEC. Furthermore, the 
present ansatz is restricted to radially symmetric sys- 
tems. To calculate excitations of cylindrically symmetric 
systems with arbitrary angular momenta, which would be 
necessary, e.g., for condensates with dipole-dipole long- 
range interactions, an extension of the ansatz is required. 
For dipolar condensates such an ansatz would be of inter- 
est, since the dipolar interaction offers the new possibility 
for a condensate to collapse with m = 2, 3, . . . symmetry, 
the so-called angular collapse [5]. 
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Appendix: Integrals for the variational ansatz with 
coupled Gaussians and spherical harmonics 

We give the integrals necessary for setting up the equa- 
tions of motion resulting from the time-dependent vari- 
ational principle for the new variational ansatz Eq. (19 1. 
We need the matrix and vector 



hi = 



dip 
dzi 
dip 
dzi 



dip 
Hip 



(A. la) 
(A. lb) 



where the mean-field Hamiltonian H consists of four 
parts 



H = T + V,, 



= -A -f 7^7-2 +8na\ip{r) 



2 / d3/M(^ 



(A.2) 
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To calculate the integrals, we write the ansatz (19) in 
a slightly different form 



For the elements of the K matrix we then obtain, with 
the definitions A^' = + (Al.)* and 7*^' = 7*^ + (7')* 



N 



(A.3) 



k=l L: 



where all cJqq = 1 have to be treated as constants, and 
not as variational parameters. 



dd] 



dtp 



1 rGi+3/2) 



dip 
' dip 



dd] 



dip 
dA^ 

dip ' 



= ~-d 



1^, T{h + 5/2l u 



(A.6) 
(A.7) 



1,, r(/, + 3/2)„,,., ^^^^^ 



(^/c/)i2+3/2 



Integrals of the K matrix 

For the elements of the K matrix, one needs the in- 
tegrals over two spherical harmonics, which because of 
their orthogonality are given by Kronecker deltas, and 
the integrals over the radial coordinate, which are all of 
the form 



dr r exp (—Ar 



(A.4) 



dip 
dM. 



dip 



dip \ 
dA^J 



dip 



dip_\ 

Qryk I 



(^H)/i+7/2 



iVw' Yd'^ r(!i+7/2).-7- 

ll .7711 

(A.9) 

9 ^ i^lrmiY^hmi 



ll ,7711 



rai + 5/2) 



(A.IO) 

I rai+3/2) 



— ^ ^ (^iimi)*^fi?7ii 



^1 I'T^ll 



(^fci)/i+3/2 



(A.ll) 



With the substitution r ^ t = Ar'^, one can use the 
Gamma function [25 to write 



Ir = ^A~'^^+^^/^T[{l + l)/2]. (A.5) 



Integrals of the kinetic term 

For the calculation of the integrals of the kinetic term, 
one lets the Laplacian act on the variational ansatz. The 
integrals of the resulting terms can then be evaluated 
using Eq. (A.5 I, which leads to 



dip 



I dip 



dip 



Tip 



N 

k 

2 / ^ "■l2m2 
fc=l 



(^fci)/2+3/2 



A{A';Y 



(^fc;)(2+5/2 



9 X! X! ('^Lmi)*^hnu 



k—1 ll ,7711 

TV 



(^fci)Zi+5/2 



(^fei);i+7/2 



k — 1 ll ,7711 



(^i;i)Zi+3/2 



(^fei)ii+5/2 



(A.12) 
(A.13) 
(A. 14) 



Integrals of the trapping potential 

The integrals for the trapping potential are straight- 
forward: 



dip 



dd\ 



l2rao 



N 



VextV' 



1 , r(;, + 5/2) . 

2 'r Z_^"-l2m2 ^J^kiy2+5/2 



k=l 



(A.15) 



\dAi^ 



1 

K=xt^ ) = -27'E E K^^rdU, 

k — 1 ll ,mi 

r{h + 7/2) 



(A.16) 
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' dip \ 1 Integrals of the scattering term 

fe-iii,mi rj,^ write down the integrals of the scattering term, 

y, r(/i + 5/2) ^_^ki introduce the new abbreviations A^,^'"'' = A"^ + A!f , 

^jjfci _ _|_ ^fci^ a,nd for the integral over four spherical 
harmonics the notation 

7^'*^ TOi; /2, m^M, ms-Ji, TO4) 

df] Yi,^, (e, m.m, {9, m^m, {0, m.n., (0, 

(A.18) 

where dfi = d0 d^ sin 9 is the differential solid angle ele- 
ment of the angular coordinates. Using again Eq. (A. 5 1, 
we obtain for the integrals 



dd] 

\. '■2 



dtp 



N 



Vih\-A7rn V V V V Yd' d'' Wi+h + h + k + i)/'2] 



— 1 /i ,mi I3 ,7713 I4 ,m4 

X In^ih, -1^2; k, -m^, I3, mg; /i, mi) e" 

N 



(A.19) 



Ana V V V V V (-l)"^+'"^(dl „, Tid^ yd^ ,„ ,„ n{h+h+h + h + 5)/2] 

(A.20) 



i,j.k—l li.mi l2;fi'i'2 13-^3 /4,m4 

X llf\l2,-m2;l4:,-m4;l3,m3;li,mi)e~'^ ' , 

N 



A \^ \^ \^ \^ \^ f T\m2+mif,l \* , jj \* n ,k ^iih + h + h + U + 3)/2] 



^i,mi ^21^2 ^3:''^3 l4,rii4 



X l'i\l2,-m2]U,-mA;h,m:i;li,mi)e ' 



(A.21) 



An analytical expression for l'^^^ is found by noting that 
the product of two spherical harmonics can be expressed 
by a series of spherical harmonics 



1=0 m=-l 



(A.22) 



where the coefficients CV^T^T^ can be written in terms 



of Wigner 3j symbols [31] 



{2h + l){2l2 + l){2l + 1) 



An 



k h I 



I 



0/ \mi 1712 



(A.23) 



Applying this expansion twice in the integral Eq. (A.18 1 
we obtain 

llf\li,mi; I2, m2] h, ms; I4, rrii) 



00 I 



=E E (-irQ-rr^' 



-mm^ 7714 
I Is I4 ' 



(A.24) 



/=0 m=-l 



The infinite sum can be cut off, since a Wigner 3j symbol 
is zero, if the triangle inequalities |/i — /2I < I < h + h 
or 1/3 — Z4I < I < I3 + h are not fulfilled, and Zi, . . . , Z4 
cannot be greater than the largest angular momentum 
included in the variational ansatz. 



Integrals of the monopolar term 



The integrals for the monopolar term read 
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9-0 
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N 



with the definition 
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(A.25) 
(A.26) 
(A.27) 



dnjdrj dfl' J dr' 
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{r'Y 



\r — r' 



(A.28) 



To calculate this integral, the monopolar interaction po- 
tential 1/ |r — r'l is expanded in terms of multipoles [25] 



^ 21 + 1^ 

1=0 m=-i > 



X 



(A.29) 



The integral /ni,p then separates into two integrals over 
the angular coordinates which can be expressed 



with the coefficients Cp^^J^" from Eq. (A.23), and one 
integral over the radial coordinates r, r', which is of Gaus- 



r 



sian type and can be solved analytically. For 
obtain 
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with the individual integrals 
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(A.30) 



(A.31) 
(A.32) 



1 



[{h+l2-l+p)/2]\ 



4 (yifei)(ii+i2-/+P+2)/2(-^j,jfe')(;3+i4+/+3)/2 ^ al \Ap'" 



2 ^ / Akl 

X rrl^) r[a3 + ?4 + ; + 2a + 3)/2]. (A.33) 



r 



The infinite sum in Eq. (A.30 1 can be cut off again due for the variational parameters for the ansatz with coupled 



to the properties of the Wigner 3j symbols. Thus all Gaussians and spherical harmonics have been calculated 
integrals necessary for setting up the equations of motion analytically. 
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